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SUMMARY: 

The implementation of traditional CFD algorithms in practical propulsion related flowfields often leads to dramatic 
reductions in efficiency and/or robustness. The present research is directed at understanding the reasons for this 
deterioration and finding methods to circumvent it. Work to date has focussed on low Mach number regions, viscous 
dominated regions and high grid aspect ratios. Time derivative preconditioning, improved definition of the local time 
stepping, and appropriate application of boundary conditions are employed to decrease the required time to obtain 
a solution, while maintaining accuracy. A number of cases having features typical of rocket engine flowfields are 
computed to demonstrate the improvement over conventional methods. These cases include laminar and turbulent 
high Reynolds number flat plate boundary layers, flow over a backward-facing step, a diffusion flame, and wall 
heat-flux calculations in a turbulent converging-diverging nozzle. Results from these cases show convergence that is 
virtually independent of the local Mach number and the grid aspect ratio, which translates to a convergence speed-up 
of up to several orders of magnitude over conventional algorithms. Current emphasis is in extending these results to 
three-dimensional flows with highly stretched grids. 

TECHNICAL DISCUSSION: 

Propulsion flowfields contain a rich variety of physical processes. The combustion process in a rocket engine takes 
place at low flow speeds and is accompanied by very high rates of energy release and species generation. Although 
the Reynolds numbers are high, mixing and diffusion are critical to achieving combustion, and recirculation regions 
are often present near the injector face. Following combustion, the flow accelerates through the nozzle to supersonic 
speeds, where the Reynolds numbers are typically very high and there is little heat release, but recombination effects are 
important for performance predictions. This wide variety of physical processes and conditions presents a formidable 
challenge for CFD codes. 

In order to design an efficient and reliable CFD code to solve these flowfields, one needs to address the issues 
which lead to a deterioration in code performance. Convergence rates and code robustness often suffer in practical 
problems due to any number of factors. For example, the initial conditions may be poorly specified, the grid resolution 
may be inadequate, very high grid stretching will generally be present, or there may be flow regimes present for which 
the algorithm was not designed to converge efficiently. An example of this latter issue occurs for traditional density- 
based algorithms in the presence of very low Mach numbers, as occurs in the combustion chamber. Density-based 
algorithms were designed for transonic and supersonic Mach number flows, but typically have difficulty computing low 
Mach number flows. The convergence rates of these algorithms decrease dramatically because of the wide disparity 

• Graduate Research Assistant 
t Research Associate 

* Distinguished Alumni Professor 


130 



between fluid particle and acoustic wave propagation speeds. Various researchers, including the present authors, have 
negated these detrimental effects through the use of time-derivative preconditioning [1-6] so that such schemes are 
now routinely used at low speeds. A related issue at low Reynolds numbers occurs where the disparity between wave 
propagation speeds and viscous diffusion rates controls the convergence. This has also been effectively addressed 
through time-derivative preconditioning [2,3,6]. 

The details of the preconditioning techniques used herein can be found in References [2] ,[3] and [7]. The basic 
technique is to premultiply the time derivative by an appropriate preconditioning matrix. For flows dominated by 
inviscid effects, the preconditioning matrix is designed to alter the acoustic speed of the governing system so that 
each of the inviscid time scales is of the same order. This allows a larger time step than the unaltered system. 
The preconditioned density-based system closely resembles a ‘pressure-based’ system at low Mach numbers, while 
reverting to the traditional density-based system at transonic and supersonic Mach numbers [8]. For viscous flows, 
there exists a time-scale disparity between the diffusive and the convective processes. At low Reynolds numbers, 
this disparity becomes strong enough to deteriorate the convergence rate. Typically the difficulty arises when the 
cell Reynolds number (e.g. Re Ax = uAx/u) becomes less than unity. Such Reynolds numbers are encountered in 
boundary layers, shear layers, recirculation regions, or anywhere the flow gradients are strong. In this case the previous 
preconditioning matrix is slightly altered so that the acoustic time scale is of the same order as the diffusive time scale 
[ 2 ]. 

It is also well known that the presence of large grid aspect ratios decreases the performance of conventional 
algorithms. Large grid aspect ratios are often required to resolve the steep gradients that occur in turbulent boundary 
layer calculations. Row calculations in long, narrow ducts also give rise to large grid aspect ratios. The decreased 
performance in the presence of large grid aspect ratios occurs because of a wide disparity in the cell crossing times in 
the different coordinate directions. The difficulties caused by the large grid aspect ratios have been addressed by the 
present authors [7], and resulted in improvements that yield convergence rates that are orders of magnitude faster than 
the conventional algorithms in highly stretched grids for two-dimensional problems. Presently, these improvements 
only apply to the 2-D ADI scheme. In order to obtain enhanced convergence, three issues must be addressed. The 
choice of the local time step should be based on the minimum CFL number 


CFLAx CFLAy 
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where A x and X y are the acoustic eigenvalues in each of the coordinate directions. The boundary conditions should be 
applied before the governing equations are approximately factored, and the viscous preconditioning should be based 


upon the maximum von Neumann number, VNN (see Ref. [7]). 


RESULTS: 

To demonstrate the efficiency and reliability of the improved algorithm, the following test cases are presented: 1) 
turbulent flat plate boundary layer, 2) H 2 /O 2 diffusion flame, 3) 2-D backward facing step, 4) turbulent rocket nozzle 
heat flux computation. Additional cases have been calculated; however, they are not presented here because of space 
limitations. 

The first case is a Mach 0.1 turbulent flat plate boundary layer computation. The Reynolds number for this case 
is 8 x 10 5 , and the maximum grid aspect ratio (located adjacent to the wall) is 8000. The k-e turbulence model was 
used for closure. Figure 1 displays the convergence history for the enhanced ADI algorithm (preconditioning + grid 
aspect ratio enhancements), and for the standard ADI algorithm. The enhanced algorithm converges in under 2000 
steps, whereas the convergence for the standard algorithm flattens after a three order of magnitude drop in the residual 
and would require over 50,000 steps to converge to machine accuracy. 
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The H 2 /O 2 diffusion flame consists of a central core flow of pure oxygen (Reynolds number of 200). The outer 
gas flow is pure hydrogen. The grid size used is 71 x 61. Both reacting and nonreacting cases were computed* and the 
convergence for both cases using the enhanced algorithm are shown in Figure 2. For the reacting case* the reaction 
rates were reduced for initial start-up and gradually increased to their full level, which is represented by the initial 
portion of the convergence curve. Once the reactions are at the full level, the convergence rate is the same as the 
nonreacting case. Figure 3 shows contours of temperature and water mass fractions for the reacting case. 

The geometry for the 2-D backward facing step corresponds to the one studied experimentally by Armaly et al. 
[9], The length and width of the duct are 10 cm and 1 cm respectively, and the step height is 0.5 cm. A uniform grid 
of 101 x 61 is used. Convergence results are shown in Figure 4 for Re = 100, 200, 300 and 400. The residuals are 
shown sequentially because each case is used as the initial condition for the subsequent case. All Reynolds number 
cases converge at about the same rate of 7 orders of magnitude in 1000 steps. Table 1 shows the predicted primary 
recirculation zone lengths, normalized by the step height. Agreement between computations and measurements is 
excellent. 

The final case shown is high Reynolds number turbulent flow through a converging-diverging nozzle with cooled 
walls. The incoming gas is at a stagnation temperature of 3500 K and the wall is maintained at 700 K. The grid size is 
141 x 121. For accuracy in heat flux predictions it is necessary to maintain a minimum of about one along the wall. 
The strong wall stretching results in grid aspect ratios that are as high as 2 x 10 6 in this region. Figure 5 shows the 
convergence for the standard algorithm and the enhanced algorithm. The enhanced algorithm converges in 2000 steps, 
whereas the standard algorithm would require over 50,000 steps to reach machine accuracy. Even though the standard 
algorithm converges about five orders of magnitude in 4000 steps. Figure 6 reveals that the wall heat flux is far from 
converged. In fact, even after 20,000 steps the wall heat flux is not quite fully converged in the throat region. This is 
because the convergence in the near wall region is much slower than the core flow due to the extremely large aspect 
ratios near the wall. Examination of the wall heat fluxes from the enhanced algorithm reveals that a fully converged 
solution is reached in just 400 iterations, indicating the flowfield in the near wall region and the core region converge 
at similar rates. 
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Figure 1: Convergence for turbulent flat plate boundary 
layer for standard and enhanced algorithm 
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Figure 2: Convergence for nonreacting and reacting 
multi-species shear layer 
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Figure 4: Convergence for backward facing step at var 
ious Reynolds numbers 


Figure 3: Contours of temperature and H^O mass frac- 
tions for the reacting shear layer 
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Figure 5: Convergence for turbulent nozzle computa- 
tions for standard and enhanced algorithms 
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Figure 6: Wall heat flux results for standard and en- 
hanced algorithms 
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